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A previously introduced multi-boson technique for the simulation of QCD with dynamical quarks is described 
and some results of first test runs on a 6 x 12 lattice with Wilson quarks and gauge group SU(2) are reported. 



1. INTRODUCTION 

The basic idea of the algorithm proposed in 
ref.[l] is to map lattice QCD to a local bosonic 
theory and to simulate the latter using standard 
techniques. The method is generally applicable, 
but attention is here restricted to the case of two 
degenerate flavours of Wilson quarks on a four- 
dimensional lattice with periodic boundary con- 
ditions in all directions. We shall first recall the 
basic steps leading to the bosonic theory and then 
present some results from test runs on small lat- 
tices with gauge group SU(2). 

2. EQUIVALENT BOSONIC THEORY 

In terms of the gauge covariant forward and 
backward lattice derivatives, V^ and V^, the 
Wilson-Dirac operator D may be written as 



D 



3 

fi=0 



{7m(V^ + V,)-V^V,}^ 



(1) 



where the lattice spacing has been set equal to 1 
for convenience. Our starting point is the effec- 
tive distribution 



Pefi[U] oc [det(D ■ 



m)fe-^^[^] 



(2) 
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of the gauge field U which one obtains after inte- 
grating over the quark fields. Sg[U] denotes the 
usual plaquette action and the bare quark mass 
m is related to the Wilson hopping parameter K 
through K = (8 -|- 2m)~^ . 

It is easy to show that the lattice Dirac opera- 
tor satisfies 

D^=j5Dj5, ||i? + m||<8 + m. (3) 

The operator 

Q = l5iD + m)/[cMiS + m)], cm > 1, (4) 

is hence hermitean and all its eigenvalues are be- 
tween — 1 and 1. In the following it will be advan- 
tageous to work with Q rather than D + m and 
so we rewrite eq.(2) in the form 

Pes[U]xdetQ^e--^^^^\ (5) 

The constant cm > 1 will be fixed later when we 
discuss the simulation of the bosonic theory. 

We now choose a sequence of polynomials P(s) 
of even degree n such that 



lim P(s) = l/s for all < s < 1 



(6) 



(an explicit example is given below). From the 
above we then infer that 



detQ^ = lim [det PiQ'-')] ^ 



(7) 



We may, furthermore, choose the polynomials so 
that their roots zi, Z2, ■ ■ ■ , z„ come in complex 
conjugate pairs with non-zero imaginary parts. 
As a consequence the root factorization may be 
written in the form 

n 

P{Q^) = constant x ]^ [(Q - ^kf + 4] , (8) 



k = l 



where 
l^k + ivk 



Zk, Vk> 0. 



(9) 



After inserting this product in eq.(7) the deter- 
minant factorizes and each factor 



det[(0-//,)2 + j.2]- 



(10) 



can be represented as a Gaussian integral over 
some auxiliary bosonic field 4>k. 

In this way one establishes the identity 



-St [[/,</.] 



(11) 



(12) 



P,n[U]= Km 1- /D[,/-]D[,/-t] 
where 

(Zi, is a normalization constant independent of 
the gauge field). 

We have thus shown that lattice QCD is a limit 
of a purely bosonic local theory. Note that the 
fields (fik are coupled to the gauge field but not 
among themselves. Their action is strictly pos- 
itive and the bosonic theory is, therefore, suit- 
able for numerical simulation. Since only a rela- 
tively small number of scalar fields can be stored 
in memory, it is however important to make sure 
that the limit n ^ cxd is reached rapidly. In other 
words, the polynomial P(s) must be chosen with 
care so as to obtain a good approximation of l/s 
already for low values of n. 

3. CHOICE OF P{s) 

The polynomial defined below approximates 
the function l/s with a uniform relative error in 
the range £ < s < 1, where e is an adjustable 
parameter. It is somewhat simpler than the poly- 
nomial proposed in ref.[l], but so far seems to do 
equally well in practice. 



Let us introduce a scaled variable u and an an- 
gle 6 through 



u = (s-e)/(l-e), 



2m- 1. 



(13) 



The Chebyshev polynomial T* (u) of degree r is 
given by 

T;(m) = cos(r6'). (14) 

For specified n and e we now define 

Pis)=[l + pT:+,{n)]/s, (15) 

where the constant p is chosen such that the 
square bracket vanishes at s = (and so is di- 
visible by s). 
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Figure 1. Relative deviation R(s) of the polyno- 
mial (15) from l/s for e = 0.05 and n = 16. 



In the range £ < s < 1 the relative fit error 



R(s) = [P(s)- l/s]s 

of this polynomial satisfies 



\R(s)\ < 2 



1 



n + l 



(16) 



(17) 



i.e. the convergence is exponential with a rate 
roughly equal to 2^/e. When s < e the polyno- 
mial continues to converge with an s-dependent 
exponential rate which approaches zero in the 
limit s ^ 0. Moreover, as can be seen from fig. 1, 
R(s) is oscillating around so that the errors may 



be expected to partly cancel when the product 
over all eigenvalues of Q'^ is taken. 
The roots of P{s) are given by 



i(l + £)cos^-iy?sm^ 



(18) 



{k = 1, 2, . . ., n). They lie on an ellipse around 
the spectral interval < s < 1 and satisfy 
Imzfc ^ (recall that n is assumed to be even). 
The transformation to the bosonic theory thus 
works out in the way explained above. 

4. NUMERICAL SIMULATION 

The bosonic theory with action S-b\[J , </)] can be 
simulated straightforwardly using local heatbath 
and over-relaxation algorithms. A few technical 
issues are addressed here and further details are 
discussed in ref.[2]. 

We first consider the updating of the "matter" 
field (/)fc . If all field variables are kept fixed except 
(/)fc at point X, the action assumes the quadratic 
form 



^k^>i>k{x.) + \<i>k{x.)^h], + constant, 



(19) 



where A-^ is an easily calculable constant positive 
matrix. The calculation of the (constant) vec- 
tor &fc requires the computation of \Q<^'k\{x) and 



\Q 9'k\{x)- To obtain the latter one also needs 
to calculate \Q<^'k\{})) at all neighbouring points y 
of X. The subroutine which applies Q to a given 
field at a given point must hence be called 10 
times. This part of the simulation thus tends to 
be rather expensive. 

There is, however, a more efficient way to pro- 
ceed. One first computes \Q<^'k\{z) at all lattice 
points z and stores the result in an auxiliary ar- 
ray. After that one passes through the lattice in 
some order and updates the field c/i^. At any given 
point X the vector 6^ is obtained easily by apply- 
ing Q to the auxiliary field. The latter must be 
corrected at z = x and all neighbouring points as 
soon as 4'k{x) has been replaced by its new value. 
This operation costs about as much as applying 
Q once at a given site. In this way the total work 



is reduced to applying Q three times per lattice 
point. 

Concerning the updating of the gauge field, we 
note that the action reduces to 



S})[U, 4>\ = ^e ii {U {x , ijl)F} -\- constant. 



(20) 



if all field variables are held fixed with the excep- 
tion of the link matrix U{x,ijl). The force F is 
a sum of two contributions, the usual one from 
the pure gauge field action and the other from 
the coupling to the "matter" fields. The widely 
known algorithms may thus be used to refresh the 
gauge field. The computation of the force how- 
ever requires an effort proportional to n, i.e. a 
gauge field update sweep tends to be as expen- 
sive as the updating of all fields 4>k. 

We finally remark that the autocorrelation 
times for heatbath updating of the fields 4>k with 
k near n/2 tend to be rather large if the parameter 
Cm in eq.(4) is set equal to 1 (as in ref.[l]). The 
probable cause for this is that the corresponding 
roots Zk are rather close to the top s = 1 of the 
spectrum of Q'^ . The latter can be lowered by 
choosing (say) cm = 1-1 and the problem then 
disappears. 

5. LOW-LYING EIGENVALUES OF Q^ 

The parameter e should be chosen such that 
most eigenvalues A of Q'^ are in the range e < 
A < 1 where -P(A) is a uniformly good approx- 
imation to 1/A (if n is large enough). To check 
whether this condition is actually met in any par- 
ticular simulation, one must compute the lowest 
eigenvalue Aq of Q^ for a representative sample of 
gauge field configurations. 

There are various methods to compute the low- 
lying eigenvalues of large sparse matrices. A par- 
ticularly safe way to obtain Aq is to minimize the 
Ritz functional 



M^) = ||0^||V||^||^ 



(21) 



in the space of all (non-zero) quark fields ip, us- 
ing a conjugate gradient algorithm [3-6]. Higher 
eigenvalues may then be computed recursively by 
minimizing the Ritz functional in the space or- 
thogonal to the approximate eigenvectors already 
found. 



It is our experience that this method works well 
in the case of the lattice Dirac operator. In partic- 
ular, the degeneracy of the eigenvalues is correctly 
obtained and the attainable relative accuracy of 
the computed eigenvalues is completely satisfac- 
tory on machines with only single precision arith- 
metic, provided one employs an accurate summa- 
tion method to compute scalar products of fields 
[7]. A rigorous bound on the final numerical accu- 
racy may actually be deduced from the stopping 
criterion of the algorithm. 

6. TEST RESULTS 

We now present some results from simulations 
of the bosonic theory on a 6^ x 12 lattice with 
SU(2) gauge group and bare parameters (3 = 2.12 
and K = 0.15. At this point the mass nip of 
the p meson is roughly equal to 1.28 (in lattice 
units), while for the pion mass m^ one obtains 
m^ = 0.93 mp. The spatial size L of the lattice 
satisfies rn^L ~ 7 so that finite volume effects are 
expected to be small. 
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Figure 2. Probability distribution of the lowest 
eigenvalue Aq of Q'^ . The full line is from sim- 
ulations of the bosonic theory with n = 60 and 
e = 0.0015, while the dashed line represents the 
result from the Hybrid Monte Carlo run. 



In our simulations of the bosonic theory a full 
update cycle involved 1 heatbath plus 1 over- 



relaxation sweep for all matter fields 4>k followed 
by 1 heatbath sweep for the gauge field [8,9]. For 
each choice of n and e we have performed at least 
128'000 and sometimes up to 192'000 such cy- 
cles using the powerful APE computers at DESY- 
IfH. The high-quality random number generator 
of refs.[10,ll] was used in these simulations and 
Cm was set equal to 1.1 in all cases. For compar- 
ision we did also run a Hybrid Monte Carlo [12] 
program for the same lattice on the CRAY-YMP 
at HLRZ. Here we generated 2610 trajectories of 
which 1610 were used for "measurement". 

The probability distribution of the lowest ei- 
genvalue Ao of Q'^ is shown in fig. 2. Note that the 
distribution should be the same for all algorithms 
which simulate the effective gauge field distribu- 
tion _Peff[?7]. The fact that there is an only small 
(and probably not significant) difference between 
the two curves in fig. 2 thus is a first indication 
that at this level of statistics one may be able to 
do with n = 60 and e = 0.0015. 

To further investigate this question we have 
made runs with e = 0.0004,0.0007,0.0015 and 
values of n ranging from 20 to 100. For e = 
0.0004 all eigenvalues A are always in the interval 
£ < A < 1, while for e = 0.0007 there are occa- 
sionally a few eigenvalues below e. It should be 
emphasized that it is not unreasonable to set e to 
a value as large as 0.0015, since -P(A) continues to 
be a good approximation to 1/A for eigenvalues A 
substantially smaller than e (cf. sect. 7). 

In fig. 3 we present our results for the plaquette 
expectation value. The data are plotted against 
the parameter 



6 = max \R(s) 

S<S<1 



(22) 



which indicates how well the polynomial P(s) ap- 
proximates l/s in the interval containing all but 
a few eigenvalues of Q^ . At fixed s and increas- 
ing n, 6 is monotonically decreasing and eventu- 
ally vanishes exponentially [eq.(17)]. In particu- 
lar, the leftmost data points in fig. 3 correspond 
to the larger values of n. 

The plot shows that there is a significant depen- 
dence of the average plaquette on n and e when 
6 is greater than 5% or so. Below this level the 
data are consistent with a constant value and also 



Table 1 

Compilation of simulation results with b < 0.05 



n 


s 


6 


(plaquette) 


(Ao) 


m^ 


nip 


60 


0.0015 


0.0177 


0.5802(6) 


0.00110(3) 


1.173(16) 


1.262(17) 


80 


0.0007 


0.0275 


0.5806(8) 


0.00107(5) 


1.186(21) 


1.257(22) 


80 


0.0015 


0.0038 


0.5801(8) 


0.00109(4) 


1.176(15) 


1.263(20) 


100 


0.0004 


0.0352 


0.5805(9) 


0.00115(3) 


1.211(20) 


1.300(20) 


100 


0.0007 


0.0095 


0.5814(10) 


0.00115(2) 


1.204(11) 


1.290(13) 


100 


0.0015 


0.0008 


0.5799(7) 


0.00114(4) 


1.208(16) 


1.290(15) 


HMC 


- 


- 


0.5796(6) 


0.00117(2) 


- 


- 
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Figure 3. Plaquette expectation value plotted as 
a function of the fit accuracy 6 [eq.(22)]. Open 
circles, full circles and crosses correspond to dif- 
ferent values of s (0.0004, 0.0007, 0.0015). The 
line is obtained by fitting the data points with 
6 < 0.05 by a constant. 



with the result obtained from the Hybrid Monte 
Carlo (cf. table 1). 

Similiar conclusions are reached when consid- 
ering the expectation value of the lowest eigen- 
value Ao of Q^ (fig. 4) and the masses of the it 
and p mesons (figs. 5,6). The latter were defined 
through 



cosh 'i 



C(T/2+l) 



(23) 



C(T/2) ' 

where C'(t) is the appropriate correlation function 
of local operators at time t and T denotes the time 



Figure 4. Same as fig. 3, but showing the expec- 
tation value of the lowest eigenvalue Aq of Q^ . 



size of the lattice (i.e. T = 12 in the present case). 
These masses are, therefore, not quite the same 
as the true masses defined through the eigenval- 
ues of the transfer matrix. To study the proper- 
ties of the new simulation algorithm the defini- 
tion (23) is, however, equally useful and perhaps 
even better suited, since extrapolation errors are 
avoided. The correlation functions C'(t) were ob- 
tained in the conventional way, i.e. by computing 
quark propagators for a representative sample of 
gauge field configurations. 

7. CONCLUDING REMARKS 

The simulations that we have performed so far 
show that the new algorithm works out in the 
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Figure 5. Results for the pion mass m^ (cf. fig. 3). 
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Figure 6. Results for the p meson mass nip 
(cf. fig. 3). 



sense that correct and accurate results can be ob- 
tained with a reasonable number of boson fields. 
Some further consolidation by studies on larger 
lattices and at smaller quark masses is of course 
needed. One would also like to better understand 
how to choose the parameters n and e in any 
given instance. For the lattice considered above 
one may take n = 60 and e = 0.0015 to obtain 
correct results at the level of statistical precision 
given in table 1. 

It is somewhat surprising that the system- 
atic effects (which are otherwise clearly visible in 



figs. 3-6) are apparently negligible for fit accura- 
cies 6 as large as a few percent. Perhaps this is 
due to a coherent cancellation of errors when tak- 
ing the product over all eigenvalues of the quark 
matrix (as already discussed in sect. 3). 

Another interesting observation is that when 6 
is smaller than about 5% the "measured" quan- 
tities did not seem to depend on the values of e 
we have considered. This suggests that the main 
effect of the quark determinant comes from the 
high-energy modes and that correlations between 
the low-lying eigenvalues in quark loops and the 
observed quantities are small. 

Before further extensive tests are performed, 
we would like to install a number of improvements 
in the present version of the simulation program. 

a. It is well-known that the condition number 
of the quark matrix can be significantly reduced 
by taking advantage of its special form with re- 
spect to the even and odd sublattices. Even-odd 
preconditioning can be incorporated in a rather 
straightforward manner in the transformation to 
the bosonic theory. As a result the same level of 
accuracy should be reached with only about half 
of the boson fields required without precondition- 
ing. 

b. As discussed in ref.[2] unexpectedly long auto- 
correlation times are observed when the bosonic 
theory is simulated following the lines of sect. 4. 
The origin of this effect is now understood and an 
idea has been put forward of how to accelerate the 
simulation algorithm [13]. 

c. Let A[U] be any observable depending on the 
gauge field U and < A > its expectation value 
with respect to the true QCD distribution Pef[[U]. 
One may then show that 



< A >=< XA >i, / < X >i, 



(24) 



where < . . . >j denotes an expectation value in 
the bosonic theory with action Si[U, (f)] and 



X[U] = det [l + R{Q^)] 



(25) 



Approximating < A > hy < A >i, (as we did 
in the test runs) thus amounts to neglecting the 
correlations between A and X. These are guar- 
anteed to be small when n is large, provided none 



of the eigenvalues A of Q'^ are substantially below 
In general we have 



x{v\ = n [1 + ^(^)] X [1 + o(^)] 



(26) 



A<£ 



and we may thus attempt to improve the con- 
vergence in the large n limit by including X in 
the "measurement" process, using eqs.(24) and 
(26). This may be particularly worthwhile in 
those instances, where there are just a few low- 
lying eigenvalues. In any case, one may in this 
way obtain further insight into the importance of 
these eigenvalues for the quantities of interest. 
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